"""
L-Shape Domain Geometry

Implements the L-shaped computational domain by removing the upper-right quadrant 
from a rectangular domain. Provides utilities for creating fluid masks, boundary 
conditions, and domain-specific operations.

The L-shape domain is defined as: Ω = [-1,1]² \\ [0,1]² (2D) or Ω = [-1,1]³ \\ [0,1]³ (3D)
"""

"""
    create_lshape_domain_impl(n::Int, ::Val{2}) -> LShapeDomain{2}

Create a 2D L-shaped domain with n interior intervals per axis.
"""
function create_lshape_domain_impl(n::Int, ::Val{2})
    # Generate coordinate arrays (including boundary nodes)
    x = range(-1.0, 1.0, length=n+2) |> collect
    y = range(-1.0, 1.0, length=n+2) |> collect
    h = x[2] - x[1]
    
    # Create coordinate meshgrids
    nx, ny = length(x), length(y)
    
    # Generate masks
    fluid = fluid_mask_2d(x, y)
    boundary = boundary_mask_2d(x, y, fluid)  
    interior = interior_mask_2d(fluid, boundary)
    
    # Domain bounds  
    bounds = ((-1.0, 1.0), (-1.0, 1.0))
    removed_bounds = ((0.0, 1.0), (0.0, 1.0))
    
    return LShapeDomain{2}(bounds, removed_bounds, fluid, boundary, interior)
end

"""
    create_lshape_domain_impl(n::Int, ::Val{3}) -> LShapeDomain{3}

Create a 3D L-shaped domain with n interior intervals per axis.
"""
function create_lshape_domain_impl(n::Int, ::Val{3})
    # Generate coordinate arrays
    x = range(-1.0, 1.0, length=n+2) |> collect
    y = range(-1.0, 1.0, length=n+2) |> collect  
    z = range(-1.0, 1.0, length=n+2) |> collect
    h = x[2] - x[1]
    
    # Generate masks
    fluid = fluid_mask_3d(x, y, z)
    boundary = boundary_mask_3d(x, y, z, fluid)
    interior = interior_mask_3d(fluid, boundary)
    
    # Domain bounds
    bounds = ((-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0))
    removed_bounds = ((0.0, 1.0), (0.0, 1.0), (-1.0, 1.0))  # Remove for all z
    
    return LShapeDomain{3}(bounds, removed_bounds, fluid, boundary, interior)
end

"""
    fluid_mask_2d(x::Vector{Float64}, y::Vector{Float64}) -> BitMatrix

Create fluid mask for 2D L-shaped domain.
Returns true for fluid nodes, false for solid (removed) nodes.

The L-shape removes the upper-right quadrant: x > 0 AND y > 0.
"""
function fluid_mask_2d(x::Vector{Float64}, y::Vector{Float64})
    nx, ny = length(x), length(y)
    fluid = BitMatrix(undef, nx, ny)
    
    @inbounds for j in 1:ny, i in 1:nx
        # Node is solid if x > 0 AND y > 0 (upper-right quadrant)  
        solid = (x[i] > 0.0) && (y[j] > 0.0)
        fluid[i, j] = !solid
    end
    
    return fluid
end

"""
    fluid_mask_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}) -> BitArray{3}

Create fluid mask for 3D L-shaped domain.
The L-shape removes the upper-right region: x > 0 AND y > 0 for all z.
"""
function fluid_mask_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64})
    nx, ny, nz = length(x), length(y), length(z)
    fluid = BitArray{3}(undef, nx, ny, nz)
    
    @inbounds for k in 1:nz, j in 1:ny, i in 1:nx
        # Node is solid if x > 0 AND y > 0 (for all z values)
        solid = (x[i] > 0.0) && (y[j] > 0.0)  
        fluid[i, j, k] = !solid
    end
    
    return fluid
end

"""
    boundary_mask_2d(x::Vector{Float64}, y::Vector{Float64}, fluid::BitMatrix) -> BitMatrix

Create boundary mask for 2D L-shaped domain.
Returns true for boundary fluid nodes where boundary conditions are applied.

Boundary nodes include:
- Exterior walls: x = ±1, y = ±1
- Cut walls: x = 0 (y ≥ 0) and y = 0 (x ≥ 0)
"""
function boundary_mask_2d(x::Vector{Float64}, y::Vector{Float64}, fluid::BitMatrix)
    nx, ny = length(x), length(y)
    boundary = BitMatrix(undef, nx, ny)
    eps = 1e-12
    
    @inbounds for j in 1:ny, i in 1:nx
        if !fluid[i, j]
            boundary[i, j] = false
            continue
        end
        
        # Check if on exterior boundaries
        on_exterior = (abs(x[i] - (-1.0)) < eps) ||  # Left wall
                      (abs(x[i] - 1.0) < eps) ||     # Right wall  
                      (abs(y[j] - (-1.0)) < eps) ||  # Bottom wall
                      (abs(y[j] - 1.0) < eps)        # Top wall
        
        # Check if on cut boundaries (internal walls)
        on_cut_vertical = (abs(x[i] - 0.0) < eps) && (y[j] >= 0.0)    # x = 0, y ≥ 0
        on_cut_horizontal = (abs(y[j] - 0.0) < eps) && (x[i] >= 0.0)  # y = 0, x ≥ 0
        
        boundary[i, j] = on_exterior || on_cut_vertical || on_cut_horizontal
    end
    
    return boundary
end

"""
    boundary_mask_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, 
                     fluid::BitArray{3}) -> BitArray{3}

Create boundary mask for 3D L-shaped domain.
"""
function boundary_mask_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, 
                         fluid::BitArray{3})
    nx, ny, nz = length(x), length(y), length(z)
    boundary = BitArray{3}(undef, nx, ny, nz)
    eps = 1e-12
    
    @inbounds for k in 1:nz, j in 1:ny, i in 1:nx
        if !fluid[i, j, k]
            boundary[i, j, k] = false
            continue
        end
        
        # Check if on any exterior boundary
        on_exterior = (abs(x[i] - (-1.0)) < eps) || (abs(x[i] - 1.0) < eps) ||
                      (abs(y[j] - (-1.0)) < eps) || (abs(y[j] - 1.0) < eps) ||
                      (abs(z[k] - (-1.0)) < eps) || (abs(z[k] - 1.0) < eps)
        
        # Check if on cut boundaries (internal walls from L-shape)
        on_cut_vertical = (abs(x[i] - 0.0) < eps) && (y[j] >= 0.0)
        on_cut_horizontal = (abs(y[j] - 0.0) < eps) && (x[i] >= 0.0)
        
        boundary[i, j, k] = on_exterior || on_cut_vertical || on_cut_horizontal
    end
    
    return boundary
end

"""
    interior_mask_2d(fluid::BitMatrix, boundary::BitMatrix) -> BitMatrix

Create interior mask for 2D domain.
Interior nodes are fluid nodes that are not on boundaries.
"""
function interior_mask_2d(fluid::BitMatrix, boundary::BitMatrix)
    return fluid .& (.~boundary)
end

"""
    interior_mask_3d(fluid::BitArray{3}, boundary::BitArray{3}) -> BitArray{3}

Create interior mask for 3D domain.
"""
function interior_mask_3d(fluid::BitArray{3}, boundary::BitArray{3})
    return fluid .& (.~boundary)
end

"""
    count_domain_nodes(domain::LShapeDomain{D}) -> NamedTuple

Count different types of nodes in the domain for diagnostics.
"""
function count_domain_nodes(domain::LShapeDomain{D}) where D
    total = length(domain.fluid_mask)
    fluid_count = sum(domain.fluid_mask)
    boundary_count = sum(domain.boundary_mask) 
    interior_count = sum(domain.interior_mask)
    solid_count = total - fluid_count
    
    return (
        total = total,
        fluid = fluid_count,
        interior = interior_count, 
        boundary = boundary_count,
        solid = solid_count,
        fluid_fraction = fluid_count / total
    )
end

"""
    is_moving_lid_node(x::Float64, y::Float64, ::Val{2}) -> Bool

Check if a 2D node is on the moving lid boundary.
The moving lid is the top boundary (y = 1) for x ≤ 0.
"""
function is_moving_lid_node(x::Float64, y::Float64, ::Val{2})
    eps = 1e-12
    return (abs(y - 1.0) < eps) && (x <= 0.0)
end

"""
    is_moving_lid_node(x::Float64, y::Float64, z::Float64, ::Val{3}) -> Bool

Check if a 3D node is on the moving lid boundary.
The moving lid is the top boundary (y = 1) for x ≤ 0 and all z.
"""
function is_moving_lid_node(x::Float64, y::Float64, z::Float64, ::Val{3})
    eps = 1e-12
    return (abs(y - 1.0) < eps) && (x <= 0.0)
end

"""
    print_domain_info(domain::LShapeDomain{D}) where D

Print diagnostic information about the domain.
"""
function print_domain_info(domain::LShapeDomain{D}) where D
    counts = count_domain_nodes(domain)
    
    println("L-Shaped Domain Information ($(D)D):")
    println("  Domain bounds: $(domain.bounds)")
    println("  Removed region: $(domain.removed_bounds)")
    println("  Total nodes: $(counts.total)")
    println("  Fluid nodes: $(counts.fluid) ($(round(100*counts.fluid_fraction, digits=1))%)")
    println("  Interior nodes: $(counts.interior)")  
    println("  Boundary nodes: $(counts.boundary)")
    println("  Solid nodes: $(counts.solid)")
end